A Numerical Renormalization Group approach to 
Non-Equilibrium Green's Functions for Quantum 
Impurity Models 

Frithjof B. Anders 

Institut fur Theoretische Physik, Universitat Bremen, P.O. Box 330 440, D-28334 
Bremen, Germany 

E-mail: anders@itp.uni-bremen.de 

Abstract. We present a method for the calculation of dynamical correlation 
functions of quantum impurity systems out of equilibrium using Wilson's numerical 
renormalization group. Our formulation is based on a complete basis set of the Wilson 
chain and embeds the recently derived algorithm for equilibrium spectral functions. 
Our method fulfills the spectral weight conserving sum-rule exactly by construction. A 
local Coulomb repulsion U > is switched on at t = 0, and the asymptotic steady-state 
spectral functions are obtained for various values of U as well as magnetic field strength 
H and temperature T. These benchmark tests show excellent agreement between the 
time-evolved and the directly calculated equilibrium NRG spectra for finite U. This 
method could be used for calculating steady-state non-equilibrium spectral functions 
at finite bias through interacting nano-devices. 
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1. Introduction 

Understanding the influence of the environment onto the non-equilibrium dynamics of 
quantum systems remains one of the challenging questions of theoretical physics. A 
finite number of quantum mechanical degrees of freedom - an orbit, a spin or a qubit 
- interacting with a infinitely large bath of non-interacting bosons or fermions with 
a continuous energy spectrum, represents a typical class of model examples for such 
systems. 

These quantum impurity models appear to be at heart of a variety of different 
physical problems. Traditionally, they were used to describe the interaction of magnetic 
impurities within a metallic host[lJ or to investigate the dissipation in quantum 
mechanics[2]. These models have contributed immensely to our understanding of 
the low temperature properties of single-electron transistors [HI H] and the tunneling 
spectroscopy of adatoms on metal surfaces. [5], [6] In addition, within the dynamical 
mean-field theory[7J [8] or its cluster extensions[9] lattice models for strongly correlated 
fermions have been mapped onto quantum impurity problems embedded in a fictitious, 
self-consistent bath. 

Many approaches to non-equilibrium are based on the Kadanoff-Baym [TO] and 
Keldysh [11] techniques. At some time t Q = a closed system characterized by a 
density operator po evolves according to the Hamiltonian TC(t). The immense difficulty 
of treating the real-time dynamics of quantum impurity systems stems from the need to 
track the full time evolution of the density operator of the entire system — environment 
plus impurity. The Kadanoff-Baym and Keldysh techniques [TTJ [10] provide an elegant 
platform for perturbative expansions of the density operator. One of the building 
blocks of such perturbative expansions are non-equilibrium Green functions. These 
non-equilibrium Green functions also contain information on the transients as well as 
the steady-state which might be reached in the long time limit for a time-independent 
Hamiltonian. In general, however, perturbative approaches are plagued by the infra- 
red divergences caused by degeneracies on the impurity, making them inadequate for 
tackling the change of ground states of quantum impurity mo dels [T2]. 

In this paper, we present a different approach for the calculation of non- 
equilibrium Green functions of quantum impurity problems. We make use of Wilson's 
numerical renormalization-group (NRG) method [T2l[T3] and its recent extension to non- 
equilibrium dynamics [Bj fT5] . Spin-spin non-equilibrium spectral functions obtained by 
a NRG calculations were investigated first by Costi about ten years ago in the context of 
the spin-boson model[16]. Here, we are interested in the evolution of fermionic spectral 
functions. We address this problem with a different approach using the complete basis 
set of the Wilson chain [T2], [13] derived in the context of the time-dependent numerical 
renormalization group [TH [15] (TD-NRG). It has already been successfully applied to 
derive sum- rule conserving equilibrium Green functions [TTJ [18]. 

We focus on a quantum impurity system characterized by the thermodynamic 
density operator p oc exp(— f3W) for times if < 0. It evolves with respect to the 
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Hamiltonian for times t' > 0. We will derive a closed analytical formula for any non- 
equilibrium Green function G(t,t') for times t, t' > given a time-independent ?v. In 
contrary to the equilibrium Green functions, [T8l [TT] transitions between different energy 
shells require a double summation over pairs of Wilson shells (m, m'). In Sec. 12.21 we 
prove that this summation can be casted into a recursion relation involving two different 
reduced-density matrices instead of the single one used in the algorithm for equilibrium 
Green functions|18[ ITT] . It can be seen analytically that only one of these two reduced- 
density matrices contributes if 7i l = W: the equilibrium algorithm [T8] is recovered. 
Therefore, the presented approach to non-equilibrium spectral functions embeds the 
equilibrium case [18], [TT] as well. 

We will heavily make use of this algorithm in another publication [19] on the 
current- volt age characteristics of interacting nano-devices. In that paper, we will derive 
a numerical renormalization group approach based on scattering states to describe 
current-carrying open quantum systems. In this formulation, the current at finite bias is 
determined by the steady-state non-equilibrium (NEQ) spectral function [20 1 12T1 1221 123] 
which depends on the density operator of the full system. At finite bias, however, the 
NEQ density operator is only known analytically for Hamiltonians which commute with 
the number operator of left and right-moving electrons. |21[ |24"] i.e. for non- interacting 
quantum impurities. This analytically known operator p Q must be evolved into the 
unknown NEQ density operator p after switching on a finite Coulomb repulsion U . 

We have used the single impurity Anderson model (SIAM) [251126] for benchmarking 
our algorithm. We have restricted ourselves to changes of local parameters of the 
quantum impurity at t = 0. Consequently, the system has evolved with respect to 
the full Hamiltonian W. For an infinitely large bath, it is expectedjlO l ITTT [24] that the 
initial p oc exp(— f37i % ) evolved into the new thermodynamic density operator of the 
fully interacting problem described by for times t —>■ oo, unless it is prohibited by 
some conservation law [24]. This is the basic underlying assumption of the perturbation 
theory in the Coulomb interaction £7j2TJ [281 121]. Therefore, the steady-state spectral 
function obtained from a time-evolved density operator should be equivalent to the 
spectra obtained directly by an equilibrium NRG calculation [25], [26], [18j [TT] . 

We will use this comparison between both spectra as benchmark for our algorithm 
in Sec. [3] We will demonstrate excellent agreement between these differently calculated 
spectral functions for switching on the local Coulomb repulsion U from U = to a finite 
value at various temperatures and local magnetic fields. 

2. Theory 

Interacting quantum dots, molecular junctions or other nano-devices are modelled by 
the interacting region TCi mp , a set of non-interacting reservoirs TCbath an d a coupling 
between both sub-systems Hi 

H = Himp + Hbath + Hi . (1) 
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We assume that the system is in equilibrium at times t < 0, and its properties 
are determined by the density operator p . One possible choice would be Tij = 0, 
which is usually the starting point of perturbative approaches based on the Keldysh 
formalism [II]. However, this is not required by our method. We only demand that the 
initial density operator can be cast in the form po = exp(—j37i l )/Z, where Ti 1 can be 
the initial Hamiltonian of the system in thermodynamic equilibrium for times t < 0. 

At t = 0, we suddenly switch from the Hamiltonian 7i = W to 7i = Ti.f. The 
retarded two-time Green function, 

G r AB (t,t') = -iTr \p [A{t + t'),B(t%] 6(t) 



iTr 



p (t')[A(t),B] s 0(f), (2) 



contains information on the correlated dynamics of two operators A and B, where 

p(t) =e~ m!t p»e mH (3) 
0{t) = e mH Oe- mH . (4) 

For fermionic operators the anti-commutator is used for B] s while for Bosonic 
operators [A(t),B] s represents a commutator. Eq. ([2]) indicates that we can interpret 
such a two-time Green function as evolving the density operator of the system from r = 
to the time r = t', and calculating the correlation function of B and A with respect 
to the relative time t > 0. We expect that when changes are restricted to the local 
part of the Hamiltonian, i. e. ?ii mp + 'Hi-, a steady-state or even a new thermodynamic 
equilibrium [101 CHI ED EI] is reached for times larger than the largest characteristic 
time-scale of the system. In these cases, the limit 

Poo = lim p(t') (5) 

exists. Eq (E]) becomes independent of t', and G(t,t') only depends of the relative time 
t in the steady-state limit. 



2.1. Complete Basis Set 

Wilson's numerical renormalization group (NRG) method is a very powerful tool for 
accurately calculating equilibrium properties of quantum impurity models. Originally 
developed for treating the single-channel, single-impurity Kondo Hamiltonian [30, 12J, 
this non-perturbative approach was successfully extended to the Anderson impurity 
model[25l [26] , and to the two-channel Anderson [31] and Kondo Hamiltonians [32l [33] . 
Recently, it was extended to equilibrium properties of impurity models with a bosonic 
bath [341 [35] . non-equilibrium dynamics of the spin-boson model [T5l [36] or even 
combinations of both fermionic and Bosonic baths [37]. 

At the heart of this approach is a logarithmic discretization of the continuous bath, 
controlled by the discretization parameter A > 1; the continuum limit is recovered for 
A — > 1. Using an appropriate unitary transformation. [12] the Hamiltonian is mapped 
onto a semi-infinite chain, defined by a sequence of finite-size Hamiltonians 7i m with 



Non-Equilibrium Green's Functions for Quantum Impurity Models 



5 



the impurity coupled to the open end. The iterations are terminated at a finite value of 
m = N which defines the Wilson chain of finite length N. The finite- size Hamiltonian 
Tim act only on the first m chain links of the Wilson chain. The length N also determines 
the temperature Tjy oc A~ N ^ 2 for which the spectral functions are calculated. For a 
detailed review on this method see Ref. [T3] . 

Recently, a complete basis set for such a Wilson chain of length iV has been 
identified [H], [15]. The set of eigenstates of 7i m can be formally constructed from the 
complete basis set {|a; mp , aco, ■ ■ ■ , a at)} of the NRG chain of length N where the a* 
label the configurations on each chain link i. Since 7i m does not act on the chain 
links m + 1, • - ■ , N, an eigenstate |r) is written as |r, e; m) where the "environment" 
variable e = {a m+1 , ■ ■ ■ ,Qn} encodes the N — m site labels a m+1 , • • ■ , a^. The index 
m is used in this notation to record where the chain is partitioned into a "subsystem" 
and an "environment" . After each iteration the eigenstates of 7i m states are divided in 
"discarded" and N s "kept" states. The standard NRG proceeds to next iteration m + 1 
using only the kept states. It was proven [TU [15] that the discarded states from all NRG 
iterations, i.e {\l, e; m)^} also form a complete basis set. Regarding all eigenstates of 
the final NRG iteration as discarded, one can formally write the Fock space of the iV-site 
chain in the form Tn = span{|Z, e; m)^}, and the following completeness relation holds: 



Here the summation over m starts from the first iteration m min at which a basis-set 
reduction is imposed. All traces below will be carried out with respect to this basis set. 
Hence, the evaluation of the spectral functions will not involve any truncation error. 
Note also that we made no reference to a particular Hamiltonian 7i in constructing the 
basis set {\l, e; m) dis \. 

At each iteration m, the Fock space JF^ of a Wilson chain with fixed length iV is 
partitioned by all previously discarded states 



N 




(6) 



m=m s 



l,e 



m—1 




(7) 



and all states present r at iteration m 



N 




m'=m l',e' 




(8) 



We will make extensive use of the completeness relation 



+ 



(9) 



rn 



in the following section. 



Non-Equilibrium Green's Functions for Quantum Impurity Models 



6 



2.2. Derivation of the NRG non- equilibrium Green function 

For the moment, we will consider only the first term of the commutator of the retarded 
Green function I(t',t) = Tr p(t')A(t)B . If the operator O t = A(t)B were a "local" 
operator, i.e. an operator which only acts on impurity degrees of freedom or a Wilson 
chain of length m min up to which all states are still maintained, we could use the TD- 
NRGpH [15] to calculate the time evolution of O t (t') = Tr p(t')O t . 

In general, the time evolution of a local operator O leads to an operator O t which 
acts on all chain degrees of freedom. Each operator O t can always be expanded in 
outer products of all many-body states spanning the Fock-space. Here, we will restrict 
ourselves always to a many-body Fock-space basis which is an approximate eigenbasis of 
the Wilson chain Hamiltonian. For the application of the TD-NRG, we require that the 
matrix elements of O t remain diagonal in and independent of the environment degrees 
of freedom e, e' 

(r,e;m\6 t \s,e';m} = 5 e , e ,0™(t) . (10) 

Then the operator qualifies as local operator as defined in Eqn. (21) of Ref. |15j . 
We insert the completeness relation Eq. ([9]) between A(t) and B and obtain the two 
contributions 

(r, e; m\A(t)B\s, e'; m) = (r, e; m\A(t)(l+ + l~)S|s, e'; m) 
= 2_,( r i e;m\A(t)\k, e"; m) (k, e";m\B\s, e';m) (11) 

k,e" 

m—l 

+ E E(^e;m|i(t)|/",e";m")^^(/ /, ,e // ;m"| J B| S ,e';m) . 

Restricting the operators A and B to local operators, the first term remains diagonal 
in e,e'[15]. In the second term, we again make use of Eq. (Q, but partitioning the 
Fock-space of the Wilson chain with respect to iteration m": 

(r,e;m|(l+„ + l m „)A{t)\l\e"-m") dls dis (l" , e"; m"\B(l+„ + l m „)\s,e';m) 
= (r,e;m\l+„Am",e";m") dts dis (l", e"; m"\B l+„\s,e';m) 
= ^2 ^(r,e;m\k 1 ,e 1 ;m"){k 1 ,ei;m"\A(t)\l'',e";m") dis 

ki,ei k2,e2 

*dis (I", e"; m"\B\k 2 , e 2 ; m"){k 2 , e 2 ; m"\s, e'\ m) 



It „ // v 



^ ^eMk^e^A^e^ 

fcl,ei fc2,e2 

x B™" k2 5 e2t e»(k2,e 2 ;m"\s,e';m) . (12) 



Note that 1~„ \s, e'\ m) = holds for m" < m, and the indices k\ and k 2 include all states 
present at iteration m" as seen from the definition of 1^„ in Eq. (jSj). The locality of the 
operators A and B has been used and leads to the condition e\ = e 2 . Since m" < m, we 
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can partition the environment degrees of freedom e\ into e\ = (e^e'J where e[ labels 
the Wilson chain degree of freedom starting from chain link m+ 1. We obtain only non- 
zero matrix elements (r, e; m\ki, e\\ m")(k2, e\\ m"\s, e'; m), if e = e[ = e'. Therefore, 
Eq. (fTOj) holds, and the matrix elements in Eq. (1121) are independent of e. 

Consequently, the operator O t = A(t)B qualifies as a local operator in the sense of 
the TD-NRG[m [15] for each time t, and I(t',t) is given by the fundamental equation 
of the TD-NRG, Eq. (3) in Ref. [H], 

TV trun 

(13) 

Here 0™ g (t) = (r, e; m\A(t)B\s, e; m) is independent of e, and reduced density matrix 

Psfr ( m ) 



e 



is given in the NRG basis of W. At each time t', the spectral information is encoded 
in the time evolution of 0(t). 

Inserting Eq. (TTTj) into Eq. (fl3]) yields two terms. The first contribution to I(t,t 2 ) 
remains diagonal in the iteration index m and is given by the following expression 

TV trun 

h(t',t)= E EE e^-^'A^e^-^ 



7TL — f^miji T^S k 

ym Jred/ 



x B™ s p^{m) . (15) 

The restricted sum Ylr™™ requires that at least one of those indices r, s labels a discarded 
state at iteration m. The second contribution to I(t',t) = Ii(t',t) + l2(t',t), l2(t',t), 
contains a double summation over the iteration indices m and m" 

TV trun m—1 

w,t)= E E E E eW ~ w 



m=m min r.s m"=m„ 



x E^H^WKVWV 



l",e" 

l'i 'I ™//| 



x dM (f,e";m"|5|s,e; m) 
x (s, e; m|p |r, e; m) (16) 

which prevents a simple evaluation of the matrix elements of A and B. Now, we insert 
Eq. (112p into Eq. ( flBT) and arrive at 

TV trun m—1 

^,t)= E E E E eW " w E<> e<(< " <) '^ 

m=m min r,s m"=m min k 1: k 2 I" ,e" 

x E^( r ' e ' m l^i> e 2! m")(s, e; m|/3 |r, e; m)(fc 2 , e 2 ; m"|s, e; m) . 

e,e 2 

(17) 
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The summation and T^" 1 ^ 1 implies that m" < m. Therefore, the 

summation can be arranged to 

N—l trun 

w,t)= E EEE^-wC 

m"=m min I" ki k2 

xplt d kl (m",t'), (18) 

where the indices k\, k 2 run over all eigenstates of H. m " present at iteration m", but the 
index I" remains restricted to the discarded states. In the last step, we have defined a 
second reduced density matrix p klk2 (m" ,t') as 

TV trun 

Pk 2 M( m "i t ') = E ^2^2( r i e ; m \ k i,ei,m") 

m=m"+l r,s e,ei 

x (k 2 , e\\ m"\s, e; m)(s, e; m\p \r, e; m) 

xe ^-EDt' . (19) 

Partitioning the environment variable e± into e\ = (a m » + i, ■ • • , a m , e'), the relation 

N trun 
m=m"+l r,a 

x (r; m|fci, {ctj}; m") (k 2 , {ctj}; m"|s; m) 



(20) 

is obtained. Here, we explicitly made use of the fact that the matrix elements 
(k 2 , ei; m"\s, e; m) are diagonal in e' and e and independent of e. The summation over 
e only enters the definition of p r s e f(m). 

Eq. (1201) connects p k e 2 d kl (rn,t') to all reduced density operators p r s ed {m') from the 
later iterations m' > m. If p k e 2 d kl (m + l,t') is given, p k2jkl (m,t') obeys the following 
recursion relation 

trun 

p r ktkA m ^)= E E ( k 2,am+i,m\s;m + 1) 

-im + l nm + 1 \ ,/ 



x [pZ d (m + ( r] m + l\k x ,a m+1 ;m) 

trun 

+ E E (* 2 ' " m+i; m i fc/; m+l ) 



-recZ 



x p r k e , k ,,(m + l,t')(k ;m+ l\k u a m +i, m) . (21) 

which we have obtained from Eq. fTSOj) . We initialize this recursion with p r fS d k u(N, t') = 0. 
Defining the auxiliary matrix 

p' r>s (m + 1, = PZ d (m + iy ^-E^)t' + ~ fe , kn{m + 1; ^ ( (22) 

the recursion relation ( |2T|) has the same structure as Eq. (40) of Ref . [15] . 

Note that the overlap matrix elements {k 2 ,a m+ i,m\k';m + 1) are identical to the 
matrix elements ^y a k2 as defined in Eq. (2) of Ref. pT]. Matrix elements of this type 
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(r; m\ki, m") can be evaluated directly using a product of m — m" such A-matrices 

HZJ. 

At each recursion step p' rs {m + involves two terms which contribute matrix 
elements to different sectors of p' . By construction, Pk',k"{ m + has only non-zero 
matrix elements for k' and k" being both retained states of the NRG iteration m + 1. 

The restricted sum over r and s projects out the other sectors of the matrix 
Prs(t') = PTri 171 + 1)*') + Ps,r{ m + f° r which at least one of the indices s,r labels 
a discarded state. Instead of a single reduced density matrix, we need to keep tract of 
two matrices at each iteration, namely p r r e f(m) and p r r e f(m,t f ). 

Then, the two contributions to I(t',t) read 

N trun 

I(t',t)= EE eW'^eW-W 



m=m min r,s k 

x B% s p^(m) 

N—l trun 
m=m mi „ I' ki k 2 

• (23) 
This formally requires only a single summation over m: the second summation over 
m' has been absorbed into the definition of p^f^ t') . Note that the index /' labels 
all discarded states at iteration m. Obviously, the same type of calculation must also 
be performed for the second term of the commutator in Eq. (T2]) in order to obtain all 
contributions for the Green function. Fourier transformation of Eq. ( 1231) with respect 
to t yields the spectral information of interest. 

It has to be emphasized that only energetic approximations have been made. The 
NRG truncation influences the partitioning of the states, but the completeness of the 
basis is always guaranteedrJU [15]. Therefore, the spectral sum-rule remains fulfilled 
exactly for each time t' as in the equilibrium case [18]. It is straight forward to apply 
our algorithm also to the lesser and greater Green functions G < (t,t') and G > (t,t f ) as 
discussed in Ref. [171. 



2.3. Steady-state limit 

For all systems in which a time-independent steady-state density operator p^ is reached, 
Eq. (jSJ) becomes equivalent to 

1 f T 

Poo = hm — / dTp (r) . (24) 



This formulation is particularly useful for a discretized representation of an infinitely 
large system since artificial finite size oscillations are averaged out. The steady-state 
limit of the two-time Green function, 

G£,(f) = Jim GXbM 

= -iTrf/U^), (25) 
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is obtained using Eqs. (1251) and ffM§ by noting that 



lun U T dre^- E ^ = 5 Er , Es . (26) 
T ^°° 1 Jo 

In the first part of Eq. fl23|) as well as in the recursion relation (j2~Tj) . the reduced density 
matrix p r s e ^(m) contributes only energy diagonal matrix elements. In general, however, 
the reduced density matrix p™%, will not be diagonal in the NRG eigenbasis. 

We introduce the integral LA,B{t') of the Fourier transformed Green function 
G r AB (uj,t / ) with respect to t as 

L A ,B(t') = - r — ^ mG A b(", • (27) 



oo 



For operators A and B, whose anti-commutator - commutator for bosonic operators - 
remains constant, L^sit') defines a sum-rule independent of t' which is fulfilled exactly 
by our approach at any time t' due to the usage of a complete basis set. Therefore, the 
averaged sum-rule 

L AtB = lim L A)B {t') = - lim - / dr —QmG r AB (uj,r) (28) 

remain exacty fullfilled as well. An example would be the single-particle spectral 
function obtained from Eq. (|23|) by setting A = f a and B = f\. In this case 
L f rt(t') = 1. In fact, we use this criterion to check explicitly the sum-rule conservation 
and found that it remains always within machine precision with an error of 10~ 15 
independent of all parameters. 

A word is in order about the usage of the term "steady-state." We expect that 
a steady-state is always reached at long times for a time independent Hamiltonian[24j 
H,f in quantum impurity systems. In a closed but infinite quantum system, where only 
T^-imp + 'Hi has been changed, the steady-state will be identical to the thermodynamic 
equilibrium described by the density operator p = exp(— flTl^/Zf, in the sense that all 
local expectation values calculated with p^ and p will be the same. It requires that the 
limit lim^oo liniy^oo is taken such that t' V in appropriate dimensionless units. 

A steady-state rather than a thermodynamic equilibrium [2 1[ [23] will be reached 
for an open quantum system in the limit t' —>■ oo [2T| 121] at finite bias. Again, it 
requires that the limit lim^oo limy^oo is taken in the correct order. However, within a 
discretized representation of such a quantum impurity system, we can never distinguish 
between the approach to a true thermodynamic equilibrium and non-equilibrium steady- 
state for times t' — > oo. Therefore, we will always use the term "steady-state" 
throughout the paper even for situations where it can be proven that the corresponding 
continuum limit of the model approaches the thermodynamic limit for infinitely long 
times[23]. In fact, the difference between our steady-state and equilibrium spectral 
function will serve as a criterion for the quality of our approach. 
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2-4- Recovering the sum-rule conserving equilibrium NRG Green function 

Equation ( |23l) must contain all contributions to the equilibrium Green function [T7l [TH] 
as well. In equilibrium, the initial and final Hamiltonian are identical (TC = 7i l = W), 
the density operator p commutes with 7i. The overlap matrix S T;S between eigenstates 
of Ti. 1 and W, S r>s = i(s]m\r,m) f must be diagonal. Then, l\ contributes with an 
energy diagonal p r s ed (N) only on the last Wilson shell and is identical to Eq. (11) of 
Ref. PS]. For m < N, p r s f(N) has only non-zero matrix elements for r and s being 
a kept state, which are explicitly excluded by the summation restriction. Therefore, 
p r s e r( m ) contributes only once to the reduced density matrix p^^ijn, t') in the recursion 
relation Eq. ( 12TT) . namely at iteration N — 1. As a consequence, the reduced density 
matrix p^^i 771 ^') becomes time independent in equilibrium and equal to the reduced 
density matrix p r k e 2 d kl (m) , i.e. Pkf^i 171 ^') = Pk^kri 171 ) ■ The Fourier transformation of 
I 2 (t' = 0,t) with respect to t yields Eq. (16) of Ref. [18]. 

2.5. The non- equilibrium NRG algorithm 

As in the equilibrium NRG, [12] each chain length N corresponds to a temperature 
T N oc A- N / 2 . For W and H f , two simultaneous NRG runs are performed in order to 
generate the density operator po using 7i % and the eigenenergies of W for the time 
evolution. At each iteration m, we calculated the overlap matrix S r y(m) between all 
eigenstates r of and all eigenstates r' of W m \\'o\. This information, as well as the 
unitary matrices diagonalizing H l m and are stored. At the end of the NRG runs, 
the equilibrium density matrix [12], [13] [HI IS] is calculated using the last iteration of 



where Z N = J2i exp(— (5 N E^). 

We have implemented the TD-NRG algorithm [15] recursively by going backwards 
from m to in — 1. For each backward iteration, we perform the following steps: 

(i) calculate the reduced density matrix in the basis of 1-C m using using Eq. (40) in 



(iv) combine pl e f(m,t') and p r r e f(m + 1) to a single reduced density matrix 

(v) evaluate the contribution of iteration m to the excitation spectrum obtained by 
Fourier transform Eq. ( J23l) 

(vi) steps (i)-(v) are repeated until we reach the iteration m m i n at which no state was 




N 



(29) 



Ref. [15] 

(ii) calculate p' rs (m + l,t') according to Eq. f T22l ) 

(iii) calculate p^(m,t') using the recursion Eq. fT2T]) 



eliminated. 



While the selection of retained states in the NRG run for 7i l is determined by the 
density matrix [i~2]. the selection of states of Jv is guided by the notion of maximizing 
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the overlap with the eigenstates of TC 1 . Amongst different truncation schemes, which 
we have implemented, the simplest was the most effective [HI [15]. In this truncation 
scheme, we selected the lowest eigenstates of at the end of each iteration m as well. 

In Ref. [19] the current through a nano-device coupled to two leads is investigated 
as function of the finite applied bias using the algorithm for NEQ spectral function 
presented here. The device is described by a two-band model. Each band representing 
the bath continuum for either left or right-moving scattering states will be set to a 
different chemical potential fi a , a = L,R. The potential different V = fiR — [il drives 
a finite current through the nano-device. In this case, the NRG run for 7i l obtains a 
faithful many-body representation of the density operator of the non-interaction problem 
(U = 0) 

p oc e -^-*>) (30) 
where operator [2T] 

a 

replaces the usual number operator for a grand canonical ensemble in order to include 
the different potentials \i a of the scattering states. 

After each iteration for TC^, one would like to retain the states with the largest 
overlap with the eigenstates of W m . These eigenstates of TL^ are generally expected to 
be connected to the eigenstates of TC 1 of the same eigenenergy relative to the ground state 
by the Lippmann-Schwinger equation for a model with a continuous bath. In practice, 
we select those eigenstates of which have the lowest diagonal matrix elements of 
the operator TC^ — Y . Therefore, the eigenenergies E s of TL^ can be divided into two 
contributions 

E S = AE S + J2v*<- (32) 

a 

The first term AE S is of the order A~ m / 2 due to the truncation scheme, and the second 
term is defined by 

(s\Y \s) = y^/itX, = y^^a(s\N a \s) . (33) 

a a 

The question of the distribution and magnitude of the the excitation energies 
AE rs = E™ — E™ entering Eq. (f23|) arises in order to understand the redistribution 
of spectral weight at finite bias. AE rs involves eigenenergies of TC^ and is given by 

AE rs = AE r -AE S + J2 K - <) (34) 

a. 

The single-particle spectral function is obtained from Eq. ( J231) by setting A = f a and 
B — /t. Only those states r and s can contribute to the spectral function whose total 
number of particles differs by exactly one electron, i.e. 
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Substituting Eq. (137)1) into (13"I|) yields the two equivalent ways of writing the excitation 
energies 



For models with a channel conservation law, | {n r a — n^) | = 0, 1 must hold. As 
a consequence, the excitation energies AE rs are centered around the two chemical 
potentials fi a . For interacting quantum impurity models which violate channel 
conservation [211 [T9]. the differences |AiV£ s | are given arbitrary numbers. By inserting 
a finite value of {n r a — n s a ) into Eq. ( l36l) or ( J37j) . it becomes apparent that the energy 
difference AE rs will be shifted away from either chemical potential by multiples of the 
chemical potential differences V = [il — /^[SH US]- 

A word is in order concerning the the frequency resolution. In the usual equilibrium 
NRG the lowest resolvable frequency [13] coincides with the temperature T/v oc Ar N l 2 
set by the length of the Wilson chain. The non-equilibrium Green functions G(t, t') 
depends on two different times. The Fourier transformation with respect to relative time 
t remains meaningful even in the limit t' — > oo, since the steady-state density operator 
Poo exists and is well defined by Eq. ( 12^1) . However, the smallest excitation energy 
resolved might be larger than uj^ m hr N / 2 due to the difference between p™ _ArRG 
obtained via Eq. (123]) and the exact steady-state density operator for a bath continuum. 
Depending on the bias V and values of U the lower boundary for frequency resolution 
increases to tui ow ~ A -m / 2 which typically m = N— I to m = N — 3. In all cases, we 
investigated in Ref. [19], the bias V remains significantly larger that ui ow . 

3. Results 

3.1. The single impurity Anderson model 

In order to demonstrate the potential of this approach, we will present results for the 
single-particle spectral functions of the single impurity Anderson model (SIAM) for 
which the equilibrium spectral functions are well studied [3B1 ESS HOj [EU HE] and can 
serve as benchmarks. 



AE, 



'rs 



AE r - AE S + (fj, R - fj, L ) (n r R - n s R ) ± /i L 
AE r - AE S + (fi L - fi R ) (n r L - n s L ) ± /jr 



(36) 
(37) 



The Hamiltonian of the SIAMjlH [25l [26] 



n 




(38) 



ka 



■imp 



Ho + Hu 



(39) 



£(«/+ 

rr ^ 





(J 
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H v =2[^ n °~ 1 ) (40) 

consists of a single local state, which we will denote with /, with energy ef and Coulomb 
repulsion U, coupled to a bath of conduction electrons with creation operators c\ a and 
energies The local level is subject to a Zeeman splitting in an external magnetic field 
H . Note that the single-particle term of the impurity Hamiltonian 7^ mp can be written 
in two different ways, i.e. the last two lines of Eq. ( |39l) which allows for a conventional 
interaction term - last line of Eq. (1391) - or non-interaction term containing the Hartree 
contribution and a particle-hole preserving interaction term Tiu [251 12S] . To obtain a 
continuous spectral function from the set of discrete 5-functions occurring in Ga,b{z), 
the occurring 5(lu— LU n ) functions are replaced by a Gaussian broadening on a logarithmic 
mesh 

5(uj - to n ) -> = exp ^ - } (41) 




buJny^T 

where b ranges typically between 0.6 < b < 1.2[42, 

The Fourier transformation of the Green function G r A B (t, t') with respect to t obeys 
the equation of motion 

zG\ B {z, = Tr [pit') [A, B] s ] + G r [HAhB (z, t') (42) 

for any time t' and a time- independent Hamiltonian W . (Note that a time-dependent 
T~tf(t) yields the usual integral equation, and Eq. (T42|) would not hold.) 

By setting A = f a and B = /J, Bulla et al. derived a simple but exact relation 
between two Green functions and the correlation self-energy [39] 

G" Jz,t') 

*%M = U f ° n r J ; (43) 

J a J a 

which is used to express the retarded Green function as 

k 

We have calculated the Green functions G r ^ NRG \(z,t') and G r j NRG \z,t') in the 
steady-state limit t' — > oo and have obtained the physical Green function via the 
equation of motion (Jl"l"j) and fT43|) . 

As long as not otherwise stated, all energies are measured in units of T = 7iV 2 p(0), 
a constant band width[T2] of p(u) = 1/{2D)Q{D - \uj\) is used with D/Y = 20. The 
number of kept states after each NRG iteration was N s = 2000. The check the accuracy, 
we calculated the sum-rule of the raw NRG spectral function by integrating the 5-peaks 
analytically and confirmed that for arbitrary parameters and number of states the sum- 
rule for the steady-state spectral function is fulfilled within machine precision of 10 -15 . 
The algorithm itself combines the time-dependent NRG [HI EE] implementation with 
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Figure 1. (color online) Comparison of the spectral function for the six different 
values of U for the symmetric case e/ = —U/2. The steady-state spectral function, 
obtain from switching Tiu — to a finite value is plotted as straight line, while the 
direct equilibrium calculation 1 18] is given by a dashed lines of same color for the same 
parameters. The inset shows the resonance in the vicinity of the chemical potential. 
The dashed line in the inset indicated the unitary limit of \/(nT). NRG parameters: 
T/D = nV 2 p /D = 0.05, A = 2, N s = 2000, b = 0.6, T -► 0. 

the calculation of the sum-rule conserving spectral functions as discussed elaborately in 
Ref. [IS]. 

3.2. Particle-hole symmetry 

3.2.1. External magnetic field H = 0. In Fig. (TJ the steady-state spectral functions for 
a particle-hole symmetric regime are compared with the equilibrium solution obtained 
directly from the standard NRG procedure [18]. In these calculations, the Hartree term 
U/2 has been absorbed into TC 1 . At time t' = 0, the Coulomb interaction Tiu is switched 
on. An excellent agreement between the equilibrium NRG result (dashed lines) and the 
long-time limit of the time-evolved spectral functions (solid lines) is found. The non- 
interacting resonant-level spectral function centered around u = evolves continuously 
into the Green function for a SIAM with finite U. The inset in Fig. [1] shows small 
deviations between the reference equilibrium spectra for 7i = and the steady-state 
spectra obtained from the Fourier-transform of Eq. (1251) in the Kondo regime. Note 
that the exponentially small Kondo scale not accessible to perturbation theories in U 
is always accounted for correctly within the NRG and, therefore, in our algorithm by 
the crossover to the fixed-point spectrum of [25, 26J. With increasing values of U 
and fixed A, the peak height decreases from its theoretical unitary limit of l/(7iT). The 
deviations are less that 1% for U = 2 and increase to approximately 11% for U = 10. 
The correct low-energy scale [25| [13] Tk proportional to the width of the resonance at 
uj = emerges as well in the steady-state spectral functions. 

We investigated also the impact of the initial level position onto the steady-state 
spectra. A different starting point for U = could be the traditional way of writing 
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Figure 2. (color online) Comparison of the spectral function for the six different 
values of U for the symmetric case e/ = —U/2. The Hartree term U/2 is absent in 
Hamiltonian of H l , and the Coulomb interaction Hjj — Un^n^ = TL^ — TC 1 is switched 
on at t' = 0. The steady-state spectral function are plotted as solid lines, while 
the direct equilibrium calculation [TB] yields the dashed lines for the same parameters. 
The colors (color online) are identical for the same values of U. The inset shows the 
resonance in the vicinity of the chemical potential. NRG parameters: as in Fig. [1] 

of the impurity Hamiltonian Hi mp = J2a ( e / — f^O f°f a + ^ n f n { which is identical to 
( 139|) . Here, the Hartree term U/2 is not absorbed into the single-particle energy and 
the Coulomb repulsion term TCu = Un^n^ is switched on at t' = 0. 

The results for this starting point are presented in Fig. EJ The steady-state spectra 
show an increasing deviation from the correct thermodynamic equilibrium spectrum 
which remains pinned at \/n for all values of U in accordance with the density 
of state sum rule[l3l H3]. All steady-state spectra remain particle-hole symmetric, 
guarantied by Hj , and the high energy feature are well reproduced. However, we observe 
deviations from the correct Abrikosov-Suhl resonance (ASR) already for moderate values 
of U > 2T. For large values of U, the ASR is almost absent in the steady-state spectra. 

The difference can be understood in the following way. By absorbing the Hartree 
term into the initial Hamiltonian 7i l , the average impurity occupation (n/) does not 
change with time. 7i % and W will flow to the same strong-coupling fixed point for 
T —>■ 0. The excellent agreement between the equilibrium reference spectrum and the 
steady-state spectrum can be seen in Fig. [TJ 

In Fig. [2J however, we have started with a non-interacting Hamiltonian which breaks 
particle-hole symmetry: the level position is located at e/ = —U /2. For increasing values 
olU/Y > 1, it corresponds to a doubly occupied level as the starting configuration 
while the final spectra must be particle-hole symmetric for ej = —U/2. The strong- 
coupling fixed point of 7i % is characterized by an additional marginal operator which is 
proportional to the strength of the particle-hole symmetry breaking[26]. For energies 
larger than the characteristic energy scale Tk, a good agreement is found for the high 
energy parts of the spectrum which is determined mainly by the mean occupation. 
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Figure 3. (color online) Comparison (a) of the majority spin spectral function for the 
six different values of U for the symmetric case e/ = —U/2 at a fixed finite magnetic 
field H — 0.2. The color coding and NRG parameters are identical to Fig. [TJ The 
steady-state spectral functions, obtain by switching Tijj = to a finite value are 
plotted as straight lines, while the direct equilibrium calculation [18] is given by the 
dashed lines with the same color for the same "HJ , In (b) U/Y = 8 and e//T = —4 has 
been kept constant while the external magnetic field is switch on. The inset shows the 
resonance in the vicinity of the chemical potential. NRG parameters: as in Fig. [1] 

However, the low energy spectrum, which contains the information on the many-body 
resonance, deviates increasingly with increasing values of U from the reference curve. 

3.2.2. Finite external magnetic field The particle-hole symmetry, present at H — 
is broken at a finite magnetic field. In Fig. [3]^a), a comparison is shown between the 
equilibrium spectral functions (dashed lines) and p(uj,f — ► oo) obtained after switching 
on a finite value of U in a fixed and finite magnetic field of H — 0.2. The position and 
height of the many-body resonance is well reproduced. The small deviations for the 
equilibrium values increase with increasing value of U. A shift in spectral weight from 
negative to positive frequencies of the majority spectrum at large values of U indicates a 
slight underestimation of the spin-polarization for values of U > 8. Due to the total spin 
conservation of the Hamiltonian, a relaxation of the total magnetization is prohibited. 
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Figure 4. (color online) Influence of the initial value of the level position in TC 1 on the 
steady-state spectrum for a fixed value of U = 8. The initial level position e/ has been 
set to £j/r = —3,-0.2,-0.1,0,0.1,0.2. The black dashed line shows the equilibrium 
NRG spectra for the small parameters as "W . The inset shows the resonance in the 
vicinity of the chemical potential. NRG parameters: as in Fig. [TJ 
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Figure 5. (color online) Comparison of the spectral function for the three different 
values of U for the asymmetric case. The initial level position £f has been set to 
Cf = 0.235,0.21,0.175 and = —2.4 The inset shows the resonance in the vicinity of 
the chemical potential. NRG parameters: as in Fig. [TJ 

This is the source of additional small deviations |14[ 115] besides discretization errors in 
the finite-size representation of the infinitely large system. 

Alternatively, we have kept U fixed and switched on a finite magnetic field H at 
t 1 = as depicted in Fig. Mjo)- Again, the equilibrium spectra is well reproduced by 
p(u, t' — ► oo). 

3.3. Particle-hole asymmetric regime 

The influence of the initial level position on the steady-state spectra is depicted in 
Fig. H]for local particle-hole asymmetric parameters el = —2.4 and U = 8. Again, we 
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Figure 6. (color online) Comparison of the steady-state spectra (solid line) for a 
fixed value of U = 8 and ej — —2.4 evolved from U = and the thermodynamic 
equilibrium spectra (dashed line) for different values of the temperature T/Y = 
0.66,0.12,0.02, 10" 3 . The initial level position e f has been set to e)/Y = 0.175. The 
black dashed line shows the equilibrium NRG spectra for the small parameters as 
TCf . The inset shows the resonance in the vicinity of the chemical potential. NRG 
parameters: as in Fig. [TJ 

start initially with U = 0. For variation of which changes the level occupancy nf 
very moderately, the steady-state spectral function shows only marginal changes. We 
observe a significant deviation from the equilibrium NRG spectral function only for a 
large negative initial value of e^/A = —3, for which the impurity is essentially doubly 
occupied. Although the shape and position of the high-energy excitation maxima are 
well reconstructed in this case, the strongly reduced spectral weight of the low frequency 
resonance close to the chemical potential requires additional spectral weight at high 
energies, a consequence of the sum-rule conserving algorithm. 

Particle-hole asymmetric spectral functions are displayed in Fig. for three different 
values of U. Here, we have chosen the non-interaction resonant level model H, 1 such that 
the low-temperature fixed point spectra is identical to the one of Tif . 

Since the algorithm always evaluates the spectral function at a finite temperature 
defined by T N oc A~ N/2 of the last NRG iteration [El [251 1251 02] we can also track the 
temperature evolution of the spectra. For one set of parameters used in Fig. El such a 
temperature evolution of the steady-state spectra is shown in Fig. El Dashed and solid 
lines of equal color (color-online) correspond to the same temperature. Fig. [6] clearly 
demonstrates that the steady-state algorithm can be used for the temperature evolution 
of spectral functions as well. 

4. Conclusion and Outlook 

We have presented a new algorithm to calculate non-equilibrium Green functions G(t,t f ) 
for quantum-impurity models. It is derived using the complete basis set for the 
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Wilson NRG chain [Tl| [To]. Therefore, the spectral sum- rule is always fulfilled exactly, 
independent of the number N s of kept states after each NRG iteration. We have shown 
the algorithm for calculating equilibrium spectral functions [171 02] is included in our 
approach for the case of an unaltered Hamiltonian 7i l = W . 

We believe, that this algorithm will open new doors for theoretical calculations 
of non-equilibrium quantum systems. In another publication [19], we have applied our 
method to a non-equilibrium problem for which the answer is not known a priori: an 
open quantum system comprising of a quantum dot coupled to two leads whose chemical 
potential difference drives a current through this interacting junction. Only for the non- 
interacting problem (U = 0), the exact solution is known[21J. However, by switching on 
the full Coulomb repulsion Tijj at finite bias, the steady-state non-equilibrium spectral 
function evolves from this initially known solution. The steady-state currents through 
an interacting nano-device is accessible to the numerical renormalization group method 
in the strong-coupling regime at finite bias. This method has the advantage that it 
is applicable to any arbitrary coupling strength, magnetic field and temperature. In 
contrast to perturbative approaches it allows the study of the crossover from the weak- 
coupling regime at high temperatures to the strong-coupling regime at low temperatures 
and finite bias. 

In this paper, we have restricted ourselves to the relevant case of switching on a 
finite Coulomb repulsion U at if — 0. Focusing on the steady-state limit if — > oo, 
we used the well studied equilibrium spectral functions of the SIAM as benchmark for 
the steady-stated spectra obtained with our method. Since a closed quantum impurity 
system will evolve into its thermodynamic equilibrium [21], if only TCi mp -\-TCi is changed, 
the deviation between the steady-state and the equilibrium spectra serves as a measure 
for the quality of the algorithm. 

We have shown that the steady-state spectral functions agree excellently with the 
corresponding equilibrium spectra even at finite magnetic field. The absorbing of the 
Hartree term into the non-interacting part of the Hamiltonian yields the best agreement 
between the steady-state spectra and the equilibrium NRG spectra directly obtained 
from 7v . The singly peaked spectrum of the resonant level model evolves into the 
typical three peak structure of the SIAM in the Kondo regime, with the lower and high 
frequency peaks resulting from charge fluctuations and a narrow many-body Kondo 
resonance emerging close to the chemical potential whose width is proportional to the 
correct low energy scale. 
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